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Abstract 


We demonstrate the power of adaptive mesh refinement with adjoint- 
based error estimates in verification of simulations governed by the steady 
Euler equations. The flow equations are discretized using a finite vol- 
ume scheme on a Cartesian mesh with cut cells at the wall boundaries. 
The discretization error in selected simulation outputs is estimated us- 
ing the method of adjoint- weighted residuals. Practical aspects of the 
implementation are emphasized, particularly in the formulation of the 
refinement criterion and the mesh adaptation strategy. Following a thor- 
ough code verification example, we demonstrate simulation verification 
of two- and three-dimensional problems. These involve an airfoil perfor- 
mance database, a pressure signature of a body in supersonic flow and 
a launch abort with strong jet interactions. The results show reliable 
estimates and automatic control of discretization error in all simulations 
at an affordable computational cost. Moreover, the approach remains ef- 
fective even when theoretical assumptions, e.g., steady-state and solution 
smoothness, are relaxed. 
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Nomenclature 


A Face area 

C Constant used in estimating discretization error, C = 2 or | 
d Distance vector from cell centroid to face centroid 
£ Discretization error 

E Total energy per unit mass 

e Discretization error with respect to a uniformly refined mesh 
F Inviscid flux tensor 

G Numerical flux function 

J Scalar functional or output, e.g., lift coefficient 

n Outward pointing unit normal 

N Number of cells 

P Prolongation operator or matrix 

p Pressure 

Q Flow solution vector of conservative variables [p, pu , pv , pic, pE] T 
R Vector of residuals 

t Time 

U Flow solution vector of primitive variables [p, u, v, ic,p] T 
u,v,w Cartesian components of velocity 
V Volume 

Greek letters 
p Error indicator 

p Density 

^ Adjoint vector 

Superscripts 
a Adjoint 

H Data reconstruction from mesh with characteristic cell-size H 
Subscripts 

H Discretization on mesh with characteristic cell-size H 

h Discretization on mesh with characteristic cell-size h—\H 

c Adjoint correction 

L Linear interpolant 

TL Trilinear interpolant 

TQ Triquadratic interpolant 

w Wall 

oo Freestream 
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There is no doubt that mesh adaptation can lead to 
significant improvement in solution accuracy. . . . What 
remains in doubt is whether the current methods of mesh 
adaptation can be brought to a sufficient level of reliability 
and robustness for routine use as a predictive tool. 

T. Baker, 1997 [1] 


1 Introduction 

Despite significant progress in solution-adaptive mesh refinement [1, 2, 3, 
4, 5], verification of flow simulations remains largely a manual procedure 
that requires expert guidance [6, 7, 8]. Most of the time is consumed by 
crafting a computational mesh, checking the results and refining the mesh 
to assess and control discretization error a . As the number and complexity 
of simulations increase (for example consider flight vehicle performance 
databases involving 10 3 -10 5 cases), manual simulation verification be- 
comes impractical. Instead, simulations are verified for only a subset of 
cases and involve “best-practice” guidelines, which do not guarantee that 
the results comply with the expected standards of accuracy [9, 10]. 

The numerical accuracy of flow simulations depends inextricably on 
discretization error. In other words, numerical inaccuracy is a conse- 
quence of the mesh. Establishing credibility for complex simulations, 
therefore, requires that mesh generation and error estimation be inte- 
gral parts of each simulation. As the mesh and flow solution evolve, a 
systematic reduction in the discretization error is achieved through use 
of error estimates derived from the flow solution on the current mesh. 
The benefit is a straightforward quantitative assessment of convergence 
because a mesh refinement study is intrinsic to every case. 

The goal of most engineering simulations is to predict a handful of 
outputs, for example aerodynamic forces and moments. In such goal- 
oriented simulations, it is most efficient to focus on discretization error 
directly affecting the outputs of interest. For example, even the simple 
problem of predicting the span efficiency factor of an isolated wing in 
subsonic inviscid flow becomes prohibitively expensive if the mesh is 
refined to follow the tip vortex far downstream. As the influence of the 
vortex on span efficiency decreases with downstream distance, so should 
the cell refinement. Experience with error estimates that do not target 
outputs, in particular direct residual or truncation-error estimates, indeed 
shows that the adapted meshes are frequently inferior to those crafted 
by experts who understand the goals of the simulation. 

Remarkably, a relatively straightforward modification to residual er- 
ror estimates allows the prediction of error in outputs: the cell- wise 

a The subject of validation (modeling error) is not considered and in our experience, 
discretization errors dominate other aspects of simulation verification such as iterative 
convergence. 
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residuals are weighted by their influence on the output. This is the idea 
behind the method of adjoint- weighted residuals, where the weights are 
obtained from the solution of an adjoint equation. The result is not only 
an estimate of error in the outputs due to discretization (for example the 
error in lift), but also a cell- wise error indicator to guide mesh refinement. 
The adjoint- weighted approach was first developed within the framework 
of the finite element method [11, 12] and extended to finite volume meth- 
ods by Giles and Pierce [13], Barth [14], and Venditti and Darmofal [15]. 
The approach has been steadily refined to improve its accuracy and ef- 
ficiency [16, 17, 18, 19, 20, 21, 22], and has been used successfully to 
establish the credibility of goal-oriented simulations [23, 24, 25, 26, 27]. 

The routine use of adjoint-based error estimates for automatic sim- 
ulation verification, however, is predicated on robust mesh generation. 
This is because failures in mesh generation often require expert interven- 
tion to resolve. If a typical simulation requires five to ten adaptation 
cycles to attain sufficient output accuracy, then the construction of an 
aerodynamic performance database for a moderate range of operating 
conditions may invoke the mesh generator ten thousand times. Therefore, 
the mesh generator must be fast and failsafe for automatic verification to 
be viable in an engineering environment. One such robust approach is the 
embedded-boundary (cut-cell) method, where the mesh is constructed 
by embedding the geometry in a regular lattice of hexahedral (Cartesian) 
or tetrahedral elements [28, 29, 30, 31, 23, 32, 33]. 

The purpose of this paper is to demonstrate that the combination 
of adjoint-based error estimates with a Cartesian cut-cell method is a 
practical approach for automatic verification of steady, goal-oriented sim- 
ulations. The paper covers the development and implementation of a 
simulation verification framework previously described in [34, 35, 36] with 
additional details and improvements. The framework uses the approach 
of Venditti and Darmofal [15] to formulate reliable error estimates and the 
approach of Aftosmis and Berger [37] for incremental refinement of nested 
Cartesian cut-cell meshes. The framework emphasizes robustness and 
efficiency, in terms of both execution speed and memory requirements, be- 
cause both are central considerations in engineering and decision-making. 

We begin with a brief review of discretization error in Sec. 2. An 
estimate of the error in user-selected outputs is expressed in terms of 
adjoint- weighted residuals using the algebraic formulation of [15]. Sec- 
tion 3 presents the salient features of the flow solver, the formulation of 
the discrete adjoint equation and the implementation of the adjoint solver. 
Section 4 explains the details of the error estimation procedure and the 
formulation of a robust refinement criterion. The simulation verification 
framework is presented in Sec. 5, including a discussion of the adapta- 
tion mechanics, error control and practical aspects of the implementation. 
The results are organized in two parts. Section 6 presents a code veri- 
fication example that establishes the accuracy of the error estimate for 
guiding mesh refinement. Section 7 demonstrates examples of automatic 
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simulation verification on a sequence of three problems of increasing diffi- 
culty. Additional examples of simulation verification through use of this 
framework can be found in [38, 39, 40, 41, 42, 43, 44, 45]. 


2 Error Estimates 

2.1 Discretization Error 

Our goal is to compute a reliable 
approximation of a scalar output 
functional J(Q), for example lift 
or drag, derived from a flow so- 
lution Q that satisfies the flow 
equations 

R(Q) = 0 (1) 

such as the Euler or Navier- 
Stokes equations. To compute 
a discrete approximation of the 
functional Jh(Qh ), the domain 
is tessellated into N control vol- 
umes with characteristic cell-size 
i7, which we call the “working” 
mesh. The flow equations are dis- 
cretized and solved to satisfy a 
system of modified partial differei 



Figure 1. Example of a uniform mesh 
refinement study showing the defini- 
tion of discretization error, £ , and the 
error relative to the next mesh, e. 

dal equations 


Rif (Qif) = 0 (2) 

where Q h — [Qi, Q 2 , • • • , Q n] T is the discrete flow solution vector, e.g., 
an algebraic vector of cell-average values, and the discrete operator R# 
represents the residual vector. Similarly, J# represents the discrete oper- 
ator used to evaluate scalar functionals, e.g., the integration of pressure 
to obtain lift given the flow solution on the working mesh Q #. 

The error in the functional due to discretization is 

£ = \J(Q)-Jh(Qh)\ (3) 

This is illustrated in Fig. 1, which shows a notional mesh refinement study 
on a sequence of nested, uniformly refined meshes. We assume that er- 
rors not related to i7, such as distance to farfield, are negligible. If the 
discretization is consistent, then the approximation Jh(Qh) converges 
to the exact solution J(Q) as the number of cells in the computational 
domain is increased and the characteristic cell-size H shrinks. In prac- 
tice, direct computation of the discretization error is difficult because it 
involves the analytic partial differential equations, Eq. 1. 
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An alternate approach is to compute the error relative to the func- 
tional on the next nested mesh 


e = \Jh(Qh) ~ Jh(Qh)\ ( 4 ) 


as shown in Fig. 1. We refer to the next uniformly refined mesh as 
the “embedded” mesh with characteristic cell-size h. Assuming that 
the problem is smooth, the discretization error £ can be expressed as a 
geometric series in e. For example, the error expression for a functional 
with second-order convergence is 


oo . 

£ = V ie = t 

At Q 


i = 0 


(5) 


and £ = 2e for first-order functionals 13 . This trades the need for the exact 
solution J(Q) for the requirement that the functional be in the asymptotic 
range with a known convergence rate. Put another way, the starting mesh 
should be sufficiently fine. The key step becomes approximating Jh(Qh) 
without solving on the embedded mesh. 


2.2 Method of Adjoint- Weighted Residuals 


To derive a reliable approximation of the functional Jh(Qh)i consider its 
truncated Taylor series expansion about the working-mesh solution 

MQh) » Jh(Qh) + 9J qqP (Qh - Qf) (6) 

The algebraic vector Q^ denotes a reconstruction of the flow solution 
from the working mesh to the embedded mesh via a prolongation operator, 
Qh = PQ h. The term Jh(Qh) is the evaluation of the functional 
using the reconstructed flow solution on the embedded mesh, e.g., lift 
computation using the reconstructed state and finer boundary resolution. 
This is usually straightforward. The challenge is the explicit dependence 
on Qh in the inner-product term of Eq. 6. 

To eliminate Q^, expand the residual equation to obtain 

Rh(Qh) = o « iMQf ) + dR jQ Q " } (Qh ~ Qh) (7) 

Note that Eqs. 6 and 7 are approximate. The derivation can be made ex- 
act through use of the mean- value linearization; however, this only defers 
the use of similar approximations to obtain computable error estimates. 
Combining Eqs. 6 and 7 gives 


Jh(Qh) ~ Jh(Qh) 


dJhjQj!) 

dQ h 


~ dRh(Qh 

d Q h 


n -1 


H\ 


^h{Qh 


(8) 


b The sum of a geometric series 1 + r + r 2 + . . . is if \r\ < 1. The common ratio 
r is 1/4 for second-order functionals and 1/2 for first-order functionals. 
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which is independent of Q^. The adjoint equation is obtained from Eq. 8 
by defining the following intermediate product 

- dK h (Qjp 
dQ h 

where the vector ij; denotes the adjoint variables. Rewriting Eq. 8 with 
the adjoint variables 

MQh) ~ J h ( Rh(Qh ) (io) 

reveals that the adjoints weight the residual errors to form a correction 
term to approximate the functional on the embedded mesh. Substituting 
Eq. 10 into Eq. 4, the error expression (Eq. 5) becomes 

s « c\ J h (Q% ) - ^ Rft(Qf ) - Jh(Qh)\ (11) 

where the constant C — | for second-order functionals and C — 2 for 
first-order functionals. 

While Eq. 11 is independent of Q^, it does require the solution of 
the adjoint equation on the embedded mesh of Eq. 9). This is 
impractical because a solution of the large, linear adjoint system can be 
nearly as expensive as a nonlinear flow solution. Various strategies exist 
to circumvent this difficulty. The initial step involves solving the adjoint 
system on the working mesh 

dR#(Qif) 

9Qh 

This solution is then prolonged to the embedded mesh to estimate 
The estimate can be sharpened by additional (implementation specific) 
procedures, such as relaxation. To explain the salient features of our 
approach, we first introduce the governing equations and the numerical 
method, and then return to evaluation of Eq. 11 in Sec. 4. 

3 Governing Equations and Numerical Method 

3.1 Flow Equations 

We solve the three-dimensional Euler equations governing compressible 
flow of a perfect gas. For a finite region of space with volume V and 
surface area A , the integral form of the Euler equations is given by 

[ Q dF + </ F • n cL4 = 0 (13) 

dt Jv J A 

where Q = [p, pu, pv, pw, pE] T , F is the inviscid flux tensor and n is the 
outward facing unit normal vector. 


^ H = 


QJh^SIh) 

9Qh 


( 12 ) 


h = 


d_M Qff) 
dQh 


(9) 
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The Euler equations are solved with a 
finite- volume method on a regular Carte- 
sian mesh with embedded boundaries. 

The body geometry is specified by a wa- 
tertight surface triangulation. The vol- 
ume mesh consists of hexahedral cells, 
except for a layer of body-intersecting 
cells, or cut cells, which are arbitrary 
polyhedra adjacent to the boundaries, as 
illustrated in Fig. 2. The mesh is viewed 
as an unstructured collection of control 
volumes to facilitate solution-adaptive re- 
finement. 

Spatial discretization uses a cell- 
centered approach, where the control volumes V correspond to the mesh 
cells and the cell- averaged value of Q, denoted by Q#, is located at 
the centroid of each cell. The control volumes are fixed in time. The 
semi-discrete form is given by 


Figure 2. A multilevel Carte- 
sian mesh with a cut-cell 
boundary. Adjacent cells can- 
not exceed 2:1 ratio. 


V H ^- + r h (Q h ) = 0 


(14) 


where Vjj is a diagonal matrix containing the cell volumes. The residual 
in each cell i is expressed as 


R. i = £ G, 

jeVi 


A j A j 


(15) 


where j denotes the j th face of volume Vi with area A , and G represents 
the numerical flux function. 

Residual evaluation for a second- 
order accurate discretization proceeds by 
linearly reconstructing the solution to 
the face centroid. This is illustrated 
in Fig. 3, for two neighboring Cartesian 
cells Z,r sharing a common face. Prim- 


t/J 

£/r 

Ui% > x # u r 

di 

d r 


itive variables, U = [p, u, v, w,p \ 1 , are 
used for the reconstruction, and the left 
and right states are given by 


Figure 3. Reconstruction of 
face midpoint value from cell 
centroids 


U L = U z + dnfcVU z U R = U r — d r (j) r VU r 


(16) 


Here d i and d r are the distance vectors from the cell centroids to the face 
centroid, VU is the solution gradient determined via a linear least-squares 
procedure and 0 is a vector of slope limiter values used to directionally 
enforce monotonic solutions [46]. The flux value at the face centroid is 
obtained via the flux-vector splitting approach of van Leer [47] 


G(Ul,Ur) = f+(U L ) + f-(U R ) 


(17) 



At the implementation level, the assembly of the residual vector is 
accomplished by a loop over the faces of the mesh. The flux contributions 
are scattered from the face and accumulated in the cells 

R/ = R/ + G A R r = R r — G A (18) 

where the sign reflects the change in the direction of the outward-pointing 
normal. 

All boundary conditions are enforced 
weakly. At the wall, zero normal velocity 
is enforced by specifying a wall flux. This 
flux is non-zero only for the momentum 
components and uses the pressure at the 
wall centroid of each cut cell (p w ) 

G w = (0, p w , .Pw, Pwi 0) (19) 

as shown in Fig. 4. Linear reconstruc- 
tion, Eq. 16, is used to compute p w . In 
the farfield, the flux function, Eq. 17, is used to compute the flux across 
faces on the boundary. The boundary state (either Ul or Ur) is set 
via Riemann invariants and linear reconstruction, Eq. 16, is used for the 
interior state. 

Steady-state solutions are obtained with a five-stage Runge-Kutta 
scheme accelerated by local time stepping and full approximation stor- 
age multigrid. The multigrid residual restriction operator is the sum 
of the residuals of the fine mesh cells enclosed by the coarse cell, while 
the prolongation operator is direct injection. Time-to-solution is further 
reduced by parallel computing using a highly-scalable domain decom- 
position scheme. For details on mesh generation and the flow solution 
algorithm, see Aftosmis et al. [48, 28, 49, 50] and Berger et al. [51]. 

We consider two classes of functionals as primary outputs of interest. 
The first are aerodynamic performance coefficients, such as coefficients 
of lift and drag, given by 

J = h — f (n • £)(p w - Poo ) dA (20) 

Qoo ^-ref J w 

where is the freestream dynamic pressure, A re f is the reference area, 
Poo is the freestream pressure and £ is the appropriate projection for 
the coefficient of interest, e.g., £ A Voo for lift. The discretization of 
Eq. 20 uses midpoint quadrature and is summed over all cut cells. The 
second class are field functionals that can be specified anywhere in the 
computational domain. An example is a line sensor for pressure 




Figure 4. Linear reconstruc- 
tion of pressure to face cen- 
troid at the wall 
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where L is the length of the line and n is a user specified exponent (usually 
1 or 2). The evaluation of line sensors involves finding the set of cells 
intersected by the line and integrating with midpoint quadrature that 
uses linear reconstruction of pressure to the line-segment midpoint inside 
each intersected cell. 

3.2 Discrete Adjoint Equation 

The adjoint equation as derived in Eqs. 9 and 12 is in discrete form, i.e. 
the discretized residual and functional operators are linearized. This is 
a consequence of computing the error relative to the embedded mesh, 
Eq. 4, instead of the exact solution, Eq. 3. For a reliable estimate of 
error via Eq. 11, the discrete adjoint solution must converge as the mesh 
is refined 

lim 'iph ^ (22) 

h — ^0 

where ip is the solution of the analytic adjoint equation obtained by 
linearizing the flow equations, Eq. 13, and functional before discretization. 
In other words, the discretization of the flow equations and functional 
must yield an asymptotically consistent adjoint discretization 0 . 

Adjoint consistency for functionals that involve wall-boundary inte- 
grals, such as Eq. 20, is prescribed by the form of the wall flux. Referring 
to Eq. 19, the transpose of the wall-flux Jacobian is column-rank deficient 
because the flux is zero for both the continuity and energy equations. For 
example, the adjoint system at the wall is given by 


o 

o 


Ip 1 


1 

1 

o 

s 

o 


Ip 2 


dJ u 

0 h x dp v fly dp v 0 




dJ v 

o 

5 s 

« s> 

5 s 

o 


1 

1 


dJ P 


where for clarity we assumed two dimensions, first-order discretization 
(p w = p) and linearized with respect to primitive variables. This reduces 
to 

fl x lp 2 + fl'y'tp 3 — 9Jp (24) 

which is a well-known analytic adjoint boundary condition derived by 
recognizing that any variation in wall-normal velocity is zero. Moreover, 
Eq. 23 shows that the boundary functional should be a function of only 
pressure, J — J(p). There are no similar restrictions on field functionals; 
these may be a function of any flow variable. 

At the farfield boundary, linearization of the flux function, Eq. 17, in 
conjunction with the Riemann invariants yields an inconsistent adjoint 
discretization for subsonic freestream conditions [52]. Since this should 
have little impact on functional accuracy if the distance to the farfield is 

c For a definition of adjoint consistency see [52, 53]; here we use this term only in 
the sense of Eq. 22. 
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sufficiently large, we leave the residuals unchanged. For an example of a 
duality preserving formulation see [53]. To avoid pollution of the error 
estimates, we omit values from cells adjacent to the farfield boundary 
and their face neighbors (two layers of cells). 


3.3 Adjoint Solver 


The solution of the adjoint equation proceeds by introducing an unsteady 
term in Eq. 12 to obtain the following semi-discrete form 


Vff + Tt a = 0 
at 

where the adjoint residual vector is given by 

R* = - dJH T 


OQh 


OQh 


(25) 


(26) 


An important consideration for adjoint solvers is memory usage relative 
to the flow solver. In general, since the adjoint solver is required to 
run on the working mesh iL, the large matrix-vector product ^ jp 
introduces a complication. The flow Jacobian ||| is a sparse matrix that 
is constant during the adjoint solution procedure. Its non-zero entries can 
be precomputed and stored, thereby minimizing the time to compute the 
matrix- vector product at each iteration. This is an effective strategy when 
dealing with implicit flow solvers that already store the flow Jacobian [54]. 
Alternatively, when dealing with explicit solvers some or all entries of 
the flow Jacobian can be recomputed when forming the matrix- vector 
product, see for example Barth [55], Giles et al. [56], Nielsen et al. [57] 
and Mavriplis [58]. 

We adopt the explicit flow-solution method outlined in Sec. 3.1 for the 
solution of the adjoint system. The matrix-vector product is recomputed 
on-the-fly at each evaluation of the residual by reusing the face-based 
approach of the flow solver. To demonstrate, consider first-order dis- 
cretization and examine the update of the flow-solver residual vector 
from an arbitrary face of the mesh as shown in Fig. 3 and Eq. 18 


R H 


+G A <- cell / 

— G A <— cell r 


(27) 


Linearize and apply the transpose operator to obtain 



1p r 


( 28 ) 
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Therefore, the adjoint residual update when sweeping over the faces of 
the mesh in a fashion analogous to the flow solver is 


OLli 

(29) 

op — T 

R* = R* + — 

(30) 


where we include the linearization of the split fluxes of Eq. 17. Note that 
the transpose reverses the operator order of the flow solution procedure. 
For example, linearizing one of the split fluxes 


di+T _ f<9uA T 
dQi ~ \dQi) \dUi) 


( 31 ) 


shows that the flux-function linearization is evaluated before the adjoint 
of the transformation from conservative to primitive variables. 

For second-order spatial discretization, the adjoint of the reconstruc- 
tion procedure, Eq. 16, requires an additional pass over the faces of the 
mesh. Linearization of the flow gradient involves only the geometry- 
dependent least-squares weights, which are already computed and stored 
by the flow solver. Gradient linearization is omitted in cut cells with 
volume fractions less than single-precision machine epsilon to avoid spu- 
rious adjoint values similar to those observed in [59]. Furthermore, the 
linearization assumes that the limiter values, </> in Eq. 16, are indepen- 
dent of the flow solution. In other words, the limiter is treated as a 
constant. Although this is mostly a pragmatic choice, in [34] we showed 
that the impact of this simplification on the accuracy of the linearization 
is relatively small. The right-hand-side of Eq. 12 is the linearization of 
the functional, i.e. Eq. 20 or 21. This linearization involves the pressure 
reconstruction procedure of Eq. 16 and is exact except for the treatment 
of the limiter. 

Since the eigenvalues of the flow Jacobian matrix are not changed 
by the transpose operator, we expect Eq. 25 to have similar stability 
properties as the flow equations, Eq. 14. Convergence to steady-state is 
accomplished using the same five-stage Runge-Kutta time marching and 
multigrid schemes of the flow solver. Giles [60, 61] derived conditions for 
Runge-Kutta time marching schemes and multigrid that ensure the same 
asymptotic convergence rate of the flow and adjoint solvers. This duality- 
preserving algorithm is implemented almost automatically, since the flow 
solver’s residual prolongation operator is a transpose of the restriction 
operator. 

Overall, the CPU time per iteration of the adjoint solver is roughly 
equivalent to the flow solver. This is because the additional cost of re- 
evaluating the matrix-vector product at each adjoint iteration is offset by 
reusing the local time step and limiter values directly from the flow solver. 
The implementation results in only a slight increase in memory usage 
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over the flow solver due to the storage of the converged flow solution and 
its gradient. Moreover, the face-based data structures and the domain 
decomposition scheme of the flow solver are reused with only minor 
modifications, see [34] for details. 

4 Error Estimation 

With flow and adjoint solutions in hand on the working mesh i7, we 
return to Eq. 10 to estimate the functional on the embedded mesh h. 
The embedded mesh is constructed explicitly and contains about 8N 
cells d . Computation of the residual R^(Q^) involves the reconstruction 
of the flow solution on the embedded mesh from the working mesh data 
= PQ h (recall Eq. 6). The value in each embedded cell is obtained 
from its parent cell by linear reconstruction, Eq. 16. This is denoted 
by Ql = PlQ h, where Pl represents the linear prolongation operator. 
No special treatment is performed at mesh refinement boundaries and 
cut cells, where irregular stencils pollute the residual on the embedded 
mesh. Instead, we rely on the adjoint weights to attenuate the residuals 
together with filters that compensate for these numerical artifacts when 
tagging cells for refinement (described in Sec. 5.3). 

Several studies [62, 15] show that higher-order reconstruction is more 
accurate for evaluating the residual on the embedded mesh. Nevertheless, 
we find that it is difficult to match the robustness of linear reconstruc- 
tion in practical applications with shocks and other strong non-linearities. 
Moreover, the implementation is straightforward since we can reuse func- 
tions from the existing flow-solver code. 

An estimate for ^ in Eq. 10 is obtained through use of the adjoint 
solution from the working- mesh. Similar to Q^, let ipff represent a 
reconstructed adjoint solution of Eq. 12 and rewrite Eq. 10 to obtain 

MQh) ~ Jh(Qh) - (V’f ) T RVQl) - tyh - th ) T «*( Ql) (32) 

V' V' 

Adjoint Correction Remaining Error 

This manipulation yields a computable adjoint- correction term that is 
generally non-zero in finite-volume methods. It can be used directly to 
obtain a better estimate of the functional or to obtain an error estimate 
via Eq. 11, as long as the last term, the remaining error , is small. This 
is likely, because the remaining error is a higher-order term. Moreover, 
its magnitude can be controlled with adaptive mesh refinement. Note 
that we tacitly assume that all higher derivatives of the functional and 
residual equations are also small. The crux becomes finding a reliable 
estimate of iph — ^ to formulate a robust adaptation criterion. 

We use a trilinear interpolant for constructing and a triquadratic 
interpolant for approximating These interpolants are based on shape 

d Nested subdivision of a hexahedron creates eight embedded hexahedra, but at cut 
cells some children may be inside the geometry. 
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functions for “brick” elements commonly used in finite-element methods. 
The use of interpolation for ^ is a compromise between accuracy, cost, 
and factors related to implementation and maintenance. On the one 
hand, interpolation reduces the quality of the remaining-error estimate 
because it models the solution error iph — with an interpolation error 
that primarily detects solution non-linearity. On the other hand, this 
approach maximizes speed since the number of arithmetic operations on 
the embedded mesh is relatively small. Furthermore, in contrast to the 
flow reconstruction, the adjoint reconstruction is not followed by residual 
evaluation, which relaxes robustness constraints in the implementation. 

Before constructing the interpolants, the adjoint solution is linearly 
reconstructed from the centroid to the vertices of each cell on the working 
mesh (including cut cells). We adopt Eq. 16 and each vertex receives con- 
tributions from all its coincident cells. The average of all contributions 
determines the vertex solution value. Put another way, this is a data 
smoothing step. Special logic is implemented at mesh refinement bound- 
aries, where the hanging vertices of small cells need additional updates 
from their big-cell neighbors. The solution at the eight vertices of each 
working-mesh hexahedron is used to form a unique trilinear polynomial 

^tl = Co + C\X + C 2 V + C3Z + c^xy + c$xz + c Q yz + c^xyz (33) 
The triquadratic reconstruction operator is given by 

^tq = Co + C\X + c 2 y + C3Z + c^xy + c$xz + c Q yz + c^xyz + c 8 x 2 
+ c$y 2 + C 102: 2 + c\\x 2 y + ci 2 x 2 z + c\^xy 2 + c^xz 2 + c\^y 2 z 

+ cieyz + cyjx yz + cisxy z + cigxyz (34) 

To determine the 20 unknown coefficients, we use the eight solution 
values (from the trilinear case) in conjunction with the solution gradient 
at the vertices. The gradient value at a vertex is determined by the 
arithmetic average of all gradients from cells common to the vertex. The 
Barth-Jespersen limiter [63] is used to prevent oscillatory reconstruction. 
The resulting over-determined system of 32 equations is solved in a least- 
squares sense. A well-behaved triquadratic interpolant is ensured by the 
addition of safeguards. These involve monitoring solution differences 
between the triquadratic, trilinear and cell-centroid values, and using the 
lower-order values when large differences are detected. 

We split Eq. 32 into an estimate for the corrected functional 

J c = Jh( Ql) - V’tq Ql) (35) 

and a cell- wise estimate of the remaining error in each cell of the working 
mesh 

Vh = ^2 “ ^tl) T R/i(Ql)j (36) 

jeVi 
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where j denotes the j th child of parent cell and r] H = [ 771 , 772 , . . . , t]n] t . 
Note that triquadratic reconstruction is used both in the functional cor- 
rection and remaining error, Eqs. 35 and 36, which is a slight departure 
from Eq. 32. This is based on the assumption that the triquadratic in- 
terpolant is our best estimate of the embedded-mesh adjoint. Hence it is 
used not only to compute the remaining error but also to get an improved 
functional estimate. 

Substituting Eq. 35 into Eq. 11 gives a computable estimate of dis- 
cretization error on the working mesh 


To define a local quantity suitable for driving adaptive mesh refinement, 
the remaining error, Eq. 36, is localized to form an error indicator \t]\h = 
[|? 7 i|, I 772 I, . . . , \vn\} T - The sum of the error-indicator values over the cells 
of the working mesh gives a bound on the estimate of the remaining error 


Since the absolute value operator prevents cell-wise error cancellation, 
this bound is quite conservative. In fact, when dealing with difficult 
simulations containing non-smooth flow and arbitrary geometry, 77 is 
typically more conservative than the value given by Eq. 37. e We return 
to this topic in the examples of Sec. 7. 

An alternate approach involves the use of the adjoint correction term, 
|«) T iMQf )|, as an error indicator [64, 18]. However, the remaining 
error term converges at about double the rate of the correction and is 
more conservative at sonic and stagnation lines, where adjoint variables 
vanish but their derivatives do not. Moreover, since the flow and adjoint 
equations are solved on the same mesh, there is an open question regard- 
ing the error indicator maintaining consistency of the adjoint solution 
as the mesh is refined. Although \t]\h is sensitive to non-linearities in 
the adjoint solution, other implementations [15, 5] use a complementary 
remaining-error term that involves an inner product of the flow solution 
with the adjoint residual. This term makes the mesh more suitable for the 
adjoint solution, but our numerical experiments did not show significant 
benefits. 

There are several extensions of the present approach for handling 
multiple outputs. In principle, each output requires its own adjoint, 
which significantly increases simulation cost. One way to reduce the 
cost is to form a discrete error equation that is solved in conjunction 
with a modified adjoint system [16]. In this work, a simpler approach is 

e Since 77 estimates the bound on the remaining error with respect to the embedded 
mesh, Crj can be used to extrapolate toward the exact value, similar to Eq. 37. 


£*C\J c -J h (Qh)\ 


(37) 
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Figure 5. Near-body cell-wise error indicator for Cd (log 10 \v\h) of airfoil 
in subsonic flow (left) and corresponding error histogram (right) 

used, where the solution of only one adjoint system is required. Multiple 
outputs are combined using a weighted-sum formula 


where K is the number of outputs and w is an array of user-specified 
constants. In practice, outputs are frequently projections of wall pressure, 
for example aerodynamic forces and moments, for which the weighted- 
sum of axial, normal and side forces works well. 

5 Mesh Adaptation 

5.1 Equidistribution of Remaining Error 

While error estimates are critical for assessing solution quality, automatic 
error control requires additional procedures that identify regions of high 
error and modify the mesh to drive the error below desired tolerances. 
Unlike traditional error indicators, such as feature detection and estimates 
of local truncation error, \rj\n is a direct (point-wise) estimate so there is 
no ambiguity regarding the selection of an indicator variable, its relation 
to the functional error and its convergence rate. 

To introduce the concept of error equidistribution, consider an ex- 
ample error map shown in Figure 5. The values of \t)\h for Cd in the 
near-body region of a Joukowski airfoil in subsonic flow (M^ = 0.4 and 
a = 1°) are shown in the left frame. A logarithmic scale is used to 
emphasize the rapid variation of the error indicator near the airfoil, with 
highest errors at the leading and trailing edges. The corresponding error 
histogram for the entire domain is shown on the right. The horizontal 
axis contains bins of the error-indicator values. High error cells lie to the 
right. The vertical axis is the percentage of cells in each bin. 

The histogram provides insight into how well the mesh fits the sim- 
ulation. In this case, the histogram is skewed, with most of the cells 
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Figure 6. Near-body adapted mesh colored by log 10 \t]\h in Cd for airfoil 
in subsonic flow (left) and corresponding error histogram (right) 


contributing little error. This indicates that the mesh is inefficient, which 
is expected because the near-body mesh is uniform. The high-error cells 
are close to the airfoil and, in particular, near regions of high curvature. 
The highest errors (above 10 -5 ) dominate 77. This is reflected by the 
location of the mean-error value, r]/N, shown as a dashed line in Fig. 5. 
The majority of cells have error of several orders below the mean, thus the 
mean is well to the right of the peak (the mode of the histogram). Hence, 
most of the computational work associated with this mesh is unnecessary. 

The basic strategy for controlling cell-size to minimize error is to 
equidistribute the error as the adaptation advances. The principle of error 
equidistribution has been demonstrated in [1, 3], among others. The goal 
is for each cell to contribute equally toward improving the accuracy of the 
simulation. Intuitively, any departure from a uniform error distribution 
implies that the mesh points could be redistributed to obtain a lower 
average error. Historically, error equidistribution has been sought in 
global error estimates, such as truncation errors or energy norms. These 
error estimates, however, tend to significantly increase computational 
cost because they may trigger refinement of all flow features everywhere 
in the domain. In contrast, adjoint error estimates seek equidistribution 
of the functional error and consequently focus only on regions (cells) 
important for predicting the functional. In other words, the adaptation 
seeks to refine cells that make \t]\h uniform and 77 small. 

Figure 6 illustrates the main idea, where an adapted mesh is generated 
for the airfoil example of Fig. 5. The adapted mesh is shown on the left 
and its histogram on the right. The histogram is nearly symmetric with 
little spread. The value of the error indicator varies less than an order 
of magnitude in over 80% of the cells. Comparing to Fig. 5, the adapted 
mesh reduces the error variation by seven orders of magnitude. All high 
error cells have migrated to the left and are clustered close to the mode of 
the histogram. Moreover, the mean is also close to the mode, suggesting 
that little computational work is wasted on low-error cells. The grading 
of the mesh is achieved through a series of adaptations and is a direct 
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reflection of the error map in each cycle, such as the one shown in Fig. 5. 
In the limit, the ideal histogram is a delta function. In practice, the 
discrete subdivision of Cartesian cells and cell-wise non-uniform error 
convergence, as well as overall computational cost, limit the tightness of 
the histogram. 


5.2 Error Control 

Efficiency of the adaptation procedure is driven by the amount of com- 
putational work needed to transform the error histogram of the initial 
mesh to a delta function and position it sufficiently far left in the region 
of low error. This can be controlled by carefully selecting refinement 
and coarsening thresholds above and below which cells are marked for 
refinement and coarsening, respectively. Choosing refinement thresholds 
where only the highest-error cells are refined yields tight histograms but 
requires many adaptation cycles to shift the mode into a region of low 
error. Therefore, the procedure should also minimize the number of cy- 
cles and, in particular, avoid solving on similar meshes when close to the 
final mesh. 

For simulation verification, perhaps the most common approach is to 
prescribe a tolerance TOL on the remaining error or directly on £. For 
example, the goal may be to construct a mesh with 77 < 0.0001 in Cd- 
This results in thresholds proportional to T0L/7V (or T0L/7V maa; ) that 
drive the mesh toward error equidistribution. In earlier work [35, 36], we 
evaluated this approach and modified it to accommodate a“worst-errors- 
first” strategy. This reduces computational cost by avoiding the problem 
of generating too many cells early in the adaptive process, before the error 
map is accurate, and then paying for these cells on every intermediate 
mesh until the highest-error cells are finally addressed in the closing 
cycles. The threshold is set to A • T0L/7V, where A > 1 is a user-specified 
array of constants that typically decrease as the adaptation advances. 

Use of the tolerance-driven approach in practice reveals several prob- 
lems. Specifying a meaningful TOL is awkward for functionals such as 
line sensors, where there is little intuition guiding reasonable choices 
of desirable error level. In addition, A is problem dependent and the 
resulting mesh growth is hard to control in difficult simulations, where 
occasional poor convergence of the flow or adjoint solver may occur and 
cause spurious error estimates in some region of the domain. An alter- 
nate approach is to minimize the error for a given cell budget, e.g., the 
goal is to construct a mesh with one million cells that predicts Cd most 
accurately. Cells are sorted on their level of error and a threshold is 
determined such that a fraction of the highest error cells is refined. The 
threshold in each cycle can be determined from statistics of the error 
distribution [65, 66 , 37], such as the mean and standard deviation, or set 
to meet a user-specified mesh growth. 
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5.3 Implementation 

The main steps of the simulation are given by Procedure 1. 


Procedure 1: Adaptive Mesh Refinement 
Input: Surface triangulation T and refinement parameter array r 
Termination Criteria: TOL, or maximum number of cells N max , 

or number of cycles M, or maximum 
level of refinement R ma x 

Hq < — I initial M esh(T) // Generate initial mesh 
for i ^ 1 to M do 

Q h < — Flow S olve(Hi-i) 

'i/jh < — Adjoint Solve{Hi- i, Qh) 
h < — EmbedM esh(Hi-i) // Uniform refinement 
r] H < — CellwiseErrorEstimate{h , i, Qh , ^h) 

v \v\h 

Hi i — Adapt Mesh{jii \v\h, Hi—\, TOL) 
break if (77 < TOL || iV > A^ max || R > 

Rmax) 

end 

Qh < — FIowSoIvc{Hm) II Optional solve 
Result: {J H , J c ,v}i> i = 0,...,M 


While the steps are standard, we note several places where we specialized 
the procedure: 

1. The surface triangulation T is held fixed. Consequently, it is impor- 
tant that the surface triangulation is sufficiently fine to support the 
final volume mesh, particularly in regions of high surface curvature. 

2. The exit criteria of the adaptation loop are positioned such that 
only the flow solution is computed on the final mesh. This is 
optional but effective in practice, because the primary outputs of 
the simulation come from the flow solver and it is sensible to use the 
already computed error map. In terms of simulation verification, 
the procedure is conservative because we expect 77 to decrease on 
the final mesh. 

3. Except for the final adaptation, only one level of refinement is 
added per adaptation cycle. We allow more levels on the final 
cycle because the final error map should be close to asymptotic. 
Coarsening is not considered, which improves robustness without 
significantly impacting the efficiency of steady simulations. 

4. Once the cells are marked for refinement, the tags are processed to 
enhance mesh smoothness. This includes filtering out refinement 
islands and voids, buffering tags (usually by one layer of cells) and 
enforcing 2 : 1 cell ratios. 
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(a) Setup 


(b) Pressure 


(c) Density adjoint 


Figure 7. Supersonic-vortex model problem (Mi = 2.25, 77 — 1, and 
r Q = 1.382) 


The main results of Procedure 1 are the functional value, its correction 
and remaining error for the sequence of meshes generated during the 
simulation. Convergence analysis of these quantities is used to verify the 
simulation. 

6 Code Verification 

An analytic supersonic vortex provides a model problem with a known 
solution [67] to verify the performance of the error estimation framework. 
We also verify the accuracy of the adaptive mesh refinement procedure to 
establish a benchmark for automatic simulation verification. The problem 
involves isentropic flow between concentric circular arcs at supersonic 
conditions, as shown in Fig. 7(a). The exact solution is given by 


where ai — pi — 1, pi — l/y, Mi — 2.25, 77 1, and r Q — 1.382. The 

output of interest is the integral of pressure over a portion of the outer 
arc 


as sketched in Fig. 7(a). This choice permits validation of the quadratic 
adjoint interpolant described in Sec. 2, which models the adjoint solution 
on the embedded mesh. Figure 7 shows that while the flow is smooth, 
the adjoint is not. 

All simulations are initialized with the analytic solution and solved to 
steady state without limiter. Dirichlet boundary conditions at the inlet 
are prescribed from the analytic solution; the wall boundary conditions 
are specified via Eq. 19. 
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Figure 8. Error in density versus number of cells (left), and initial mesh 
and meshes after first two uniform refinements (right) 

6.1 Uniform Mesh Refinement 

A uniform mesh refinement study is performed that involves a sequence 
of six nested meshes. Figure 8 shows the initial mesh, which contains 
47 cells. The final mesh contains 33,069 cells. Note that the boundary 
discretization changes non-smoothly with refinement, e.g., on each mesh 
the cut-cell area fractions varied by at least five orders of magnitude. 
Consequently, convergence rates are obtained via linear regression of the 
finest four solutions. 

Figure 8 shows that the solution error (point- wise Li norm of p — Ph) 
is 0(h ?) when measured over the entire domain. This is consistent with 
the second-order accurate discretization and smooth analytic solution. 
The second line in Fig. 8 shows that the convergence rate is reduced to 
slightly below 0(/i 3 / 2 ) when confined to include only the boundary due 
to the irregular discretization stencils in the cut cells. 

Despite the slower convergence of the solution along the boundary, 
the first line in Fig. 9(a) shows that error in J# is 0{h 2 ). The second 
line shows the convergence rate of error in the corrected functional J c . 
As expected, its convergence parallels J#, but the error improves by 
about a half-order of magnitude. The third line in Fig. 9(a), labeled “J c 
Exact \ is obtained by solving the adjoint equation on the embedded 
mesh and using these values in Eq. 35 when computing the corrected 
functional. This is referred to as an exact correction, which we use to 
validate the triquadratic interpolant f . Except on the initial mesh, the 
quadratic correction performs as well as the exact correction. 

Figure 9(b) examines the accuracy of the corrected functional J c in 
detail. The correction should capture the dominant part of the func- 
tional’s value on the embedded mesh J \ relative to the working mesh 
value J/f, as sketched in Fig. 1 and expressed by Eq. 4. We examine 

f Recall that the cost of an embedded-mesh adjoint solution is prohibitive in practice. 
It is roughly equal to that of a flow solution on the embedded mesh. 
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(a) Error 


(b) Relative error 


Figure 9. Convergence of functional (Jh) and its correction (J c ) for 
the uniform mesh refinement study of the supersonic vortex. Error is 
measured with respect to the exact value (left) and the value on the 
embedded mesh (right). 


the accuracy of the corrected functional using both the quadratic adjoint 
interpolant and the embedded- mesh adjoint, labeled as “Exact t/V’ in 
Fig. 9(b). Both corrections are superconvergent; the exact correction 
attains almost 0(h 3 ) when measured in this relative-error metric. The 
quadratic mimics the exact correction at a slightly slower rate and with 
about a half-order offset starting from the coarsest mesh. These results 
show that the adjoint field on the embedded mesh is accurately predicted 
with the limited tri-quadratic interpolant (described in Sec. 4). 

The ability to compute relative error accurately suggests that the 
true error should also be predicted reliably as H tends to zero. Since the 
convergence rate of the functional is second-order, Eq. 37 with C = | 
estimates the true error. The first line in Fig. 10(a) shows convergence of 
the true error | J — Jh | , which is copied from Fig. 9(a) for reference. The 
second line (squares) is the error estimate £ using Eq. 37. The estimate 
is sharp; the ratio of the estimate to the true error (effectivity) is close 
to one on all but the coarsest meshes, i.e., the lines essentially over-plot. 

The vanishing gap between the estimated and true error is quantified 
by the first line in Fig. 10(b), labeled “Exact”. We use this error gap 
to test the accuracy of the remaining error term of Eq. 32. This is 
done by evaluating Eqs. 36 and 38 without localization, i.e. by omitting 
the absolute value operator in Eq. 38. The second line in Fig. 10(b), 
labeled “Estimate”, shows the value of 4/3?7 without localization 8 . The 
estimate is within half-order of magnitude of the measured gap and has a 
similar convergence rate. The agreement is excellent considering that the 
quadratic interpolant is used in both the correction and remaining error 

g The factor of 4/3 is applied to extrapolate the remaining error to its analytic value 
based on the observed second-order functional convergence. 
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(a) Error estimate 


(b) Remaining error 


Figure 10. Accuracy of the error estimate using Eq. 11 (left). Remaining 
error and bound estimate (right) for the uniform mesh refinement study 
of the supersonic- vortex problem. 


terms, and only linear prolongation of the flow solution is used when 
evaluating the embedded- mesh residuals. The third line in Fig. 10(b), 
labeled “Error Indicator” , shows the effect of localization when forming 
the error-indicator bound 77 using Eq. 38. Recall that the absolute value 
operator prevents cell- wise error cancellation. As a result, the bound 
is roughly one-order of magnitude larger than the remaining error on 
the coarsest mesh and its convergence rate is slower, indicating that the 
bound is a conservative estimate of the remaining error. 

Having established the accuracy of the error indicator, Fig. 11 shows 
histograms of the error indicator \t]\h on the initial and final meshes. 
The histograms are similar to the ones presented in Figs. 5 and 6 except 
we use a log 2 scaling of the cell-wise error. This scaling is intuitive for 
predicting how far a bin moves to the left once its cells are refined because 
the subdivision of Cartesian cells is discrete. For example, we show in 
Fig. 10(b) that the convergence rate of the error-indicator bound 77 is 
roughly 0(/i 2 ), which implies that the mean error should shift four units 
to the left after each refinement 11 . Figure 11 shows that the mean error is 
close to 2 -12 on the initial mesh and shifts to 2 -32 after five refinements. 

Figure 11(b) shows that the uniformly refined mesh is not particularly 
well suited for the simulation at hand. The error distribution is far from 
a delta function and more than a third of the 33,069 cells contribute no 
error. The adjoint field shown in Fig. 7(c) provides some insight. The 
large variations reveal the influence of point-source mass perturbations, 
which include interactions with the inner arc, on J (Eq. 41). Recall that 
J is defined over about 2/3 of the outer arc, at which point the adjoint 
variable vanishes because any perturbation downstream of this location 
cannot influence J in this hyperbolic problem. Hence, cell refinement 

h The mean shifts p + d units to the left, where p is the order of the error indicator 
and d is the spatial dimension. 
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(a) Initial mesh (b) Final mesh (cut-cell percentage magni- 

fied 6x for clarity) 


Figure 11. Histograms of the cell- wise error indicator \ tj\h for the uniform 
mesh refinement study of the supersonic vortex. Heights are normalized 
by the total number of cells. Error indicator values in cut cells are shown 
as rectangles with white fill. 


outside the functional’s zone of dependence yields many cells with zero- 
error contributions. These cells increase cost without improving accuracy 
of the simulation. 

6.2 Adaptive Mesh Refinement 

Accuracy of the adaptive procedure is demonstrated by performing five 
adaptive refinements starting from the initial mesh of the uniform study. 
The refinement threshold is set to the mean error value, except in the first 
two adaptations where it is shifted two log 2 bins to the left of the mean 
due to the coarse initial mesh. Figure 12 shows the error convergence 
of the functional and 77 , and compares them with those of the uniform 
meshes. We observe that after each refinement the functional error and 77 
nearly match the values obtained from the uniform meshes but use fewer 
cells. The final mesh contains 13,929 cells compared to the 33,069 cells 
of the finest uniformly refined mesh. 

Figure 13 shows the final mesh and its error histogram. The refine- 
ment pattern is clearly driven by features of the adjoint solution, recall 
Fig. 7(c). Essentially no mesh refinement occurs past the functional’s 
zone of dependence. An error indicator based on local-truncation errors 
would unnecessarily adapt cells in this region. The largest sensitivity 
to residual errors occurs along the inner arc, which may seem counter- 
intuitive. This is due to the nonuniform inlet Mach number and the 
number of local Mach wave reflections that can reach the functional from 
the inner arc. More importantly, the final mesh appears to strike a good 
compromise for accommodating both the flow and adjoint solutions. The 
inner product within the error indicator \ ij\h is dominated by the ad- 
joint interpolation error — ^tl), which emphasizes regions of high 

adjoint curvature that are captured in the final refinement pattern. 
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(a) Functional error 


(b) Error-indicator bound 


Figure 12. Convergence of the functional error (left) and error indicator 
r] (right) on adapted meshes of the supersonic vortex problem. 
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Figure 13. The final adapted mesh (left) and the corresponding error- 
indicator histogram (right) for the adaptive refinement of the supersonic 
vortex problem. The cut-cell percentage is inflated 3x for clarity. 
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The histogram in Fig. 13 shows that the final adapted mesh fits the 
simulation well and contrasts markedly from that in Fig. 11(b). The 
histogram is symmetric with little spread, the mean error is close to 
the mode and the adapted mesh has almost no cells with zero error. 
Moreover, the histograms of Figs. 11 and 13 show not only the error 
indicator values in every cell, but also isolate the errors in the cut cells. 
Recall that Fig. 8 shows slower convergence rate of the solution in the 
cut cells. To compensate, the error indicators in cut cells are inflated by 
a factor of 1.5 before selecting cells for refinement. This value is based 
on experience with many problems. Figure 13 shows that mode of the 
error distribution in cut cells roughly corresponds to the mode for all 
cells, indicating that the adaptation mechanics are not adversely affected 
by the non-uniform error convergence. 


7 Examples of Simulation Verification 

Three examples are presented to demonstrate the effectiveness of Proce- 
dure 1 in simulation verification. The examples progressively increase 
in complexity, from airfoil simulations to three-dimensional simulations 
of a launch abort vehicle. To varying degrees, the examples deliberately 
violate the assumptions made in the derivation of the error estimates, in 
particular assumptions of smoothness, steady-state and fineness of the 
initial mesh. These assumptions rarely apply in practical engineering sim- 
ulations and our goal is to characterize the performance of Procedure 1 
in such situations. 

7.1 Airfoil Performance Database 

The first example demonstrates 
verification of a model aerody- 
namic performance database. The 
goal is to predict the drag coeffi- 
cient, J — Cdi of the familiar 
NACA 0012 airfoil over a range 
of freestream conditions. We con- 
sider a total of 60 cases involv- 
ing six subsonic Mach numbers, 

= {0.1,0.2,0.3,0.5,0.7,0.9} 
and four supersonic Mach numbers, 

Afoo = {1.1, 1.3, 2, 3}, at six angles 
of attack, a — {0, 0.5, 1, 2, 4, 8°}. 

Figure 14 shows the near-body 
region of the initial mesh, including a closeup of the NACA 0012 airfoil 1 

lr The trailing edge of the airfoil is modified to be sharp, which is accomplished by 
modifying the last coefficient of the equation defining its thickness distribution [68] . 



Figure 14. Near-body view of 
initial mesh with inset showing 
NACA 0012 airfoil 
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Figure 15. Cd for all freestream conditions of the airfoil performance 
database example 


with a unit chord. The distance to the outer boundary is 64 chords. The 
initial mesh contains 16 x 16 cells with characteristic length H — 8 chords. 
The airfoil is intersected by just four cells. The spatial discretization 
does not involve a limiter when M ^ <0.5; cases with M ^ >0.5 use the 
Van Leer limiter. 

To construct the database, we set the maximum level of refinement, 
Rmaxi to 18 and set the adaptation threshold to the mean error. This 
is a compromise of the various database strategies we examined in [39] , 
which contrasted uniform-error and fixed-mesh databases. In general, 
specifying a uniform error tolerance, e.g., £ < 0.0001 in C^, for all 
cases is not cost-effective. This is because the cost of a constant-error 
database is frequently dominated by a few corner cases, for which a less 
accurate computation would suffice. Alternatively, specifying the number 
of adaptation cycles (M in Procedure 1) and the mesh growth for each 
cycle fixes the cost of the database. This strategy, however, restricts the 
degree of control over the variation of error across the database. The 
current approach achieves a balance by requiring that the final meshes 
for all cases contain the same smallest cell-size (H « 0.0002 chords), but 
allows the total number of cells to vary. 

Figure 15 shows the final value of Cd for all simulations in the database. 
Although there are no analytic solutions for the NACA 0012 airfoil in a 
finite computational domain, we expect Cd to approach zero in a shock- 
free two-dimensional inviscid flow. Reassuringly, Fig. 15 shows that Cd 
is essentially zero when <0.5. Moreover, while the variation in Cd 
with angle of attack is generally mild at fixed M^, there is a rapid rise in 
Cd with respect to in the transonic regime. The largest values of Cd 
are observed when M ^ = 0.9, which is due to a strong expansion over 
the aft portion of the airfoil caused by shocks at the trailing edge. As 
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(a) Subsonic cases 


(b) Transonic cases (c) Supersonic cases 


Figure 16. Convergence of error indicator 77 for subsonic M ^ = 
{0.1, 0.2, 0.3, 0.5}, transonic M ^ — {0.7, 0.9} and supersonic M ^ — 
{1.1, 1.3,2, 3} cases of the database example 


the Mach number increases into the supersonic regime, there is a gradual 
decrease in Cd as the bow shock becomes more oblique. 

While Fig. 15 demonstrates that the simulation data follows expected 
trends, it offers little quantitative evidence for simulation verification. 
Since this example involves 60 simulations, we focus the current discussion 
only on the error estimate, Eq. 37, and the error indicator, Eq. 38, which 
are central to verification and adaptation. The single-point examples in 
Subsections 7.2 and 7.3 study these quantities in greater detail. 

Error-Indicator Convergence 

Figure 16 shows convergence of the error-indicator bound 77 (Eq. 38). 
The cases are divided according to the freestream Mach number into 
groups of subsonic, transonic and supersonic flow. Overall, 77 converges 
well. This implies mesh convergence for all cases and demonstrates the 
ability of the adaptation to control the magnitude of the remaining error 
in Eq. 32. The circle symbols, plotted for an arbitrary case in each plot, 
mark adaptation cycles. The error indicator initially increases until the 
meshes reach around 1000 cells. This is typical when starting from such 
a coarse mesh (recall Fig. 14). Increasing 77 indicates that new features 
of the solution are emerging, which are not yet captured by the error 
analysis. Nevertheless, the early error maps are sufficiently accurate to 
identify critical regions of the evolving flowfield. When combined with 
the incremental /i-refinement strategy of Procedure 1, the error maps 
reliably handle coarse initial meshes. As the simulations approach i? max , 
in particular over the last four meshes of each case, 77 is decreasing linearly 
indicating that the output is asymptotic. 

Figure 16(a) shows that 77 is reduced by three orders of magnitude 
for all but two of the 24 subsonic cases. The convergence pattern is 
very similar among the cases, especially once the meshes reach 1000 cells. 
The four cases with the sharp error rise at 1000 cells are all at a = 0°. 
The case with the slowest convergence is when = 0.5 and a = 8 °, 
indicated by the dashed line, which corresponds to the onset of transonic 
flow. To reach the desired Rmax of 18, the meshes for all subsonic cases 
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left) and samples of Cd convergence on the last four meshes 


require 14 adaptations^ . 

Figure 16(b) shows convergence of 77 for the 12 transonic cases. The 
meshes, on average, are larger than those for the subsonic cases. The 
cases with the smallest values of 77 and also the smallest number of cells 
correspond to the shock-free cases, i.e., M ^ = 0.7 and a < 1°. The 
slowest convergence (largest 77 ) occurs when — 0.7 and a — 8 °; the 
next slowest is M ^ = 0.7 and a = 4°. These two cases are actually 
the slowest to converge in the entire database. We obtained similar 
results in [39] , where the largest sensitivity of aerodynamic performance 
to discretization error was also near = 0.7. The primary culprit is 
the dependence of Cd on the location of the upper surface shock. Once 
Moo reaches 0.9, the shocks migrate to the trailing edge. Figure 16(c) 
shows convergence of the 24 supersonic cases. These cases converge well 
and are tightly clustered. These meshes require 15 adaptations to reach 
Rmax of 18, which is one additional cycle when compared to the subsonic 
and transonic cases. 

Level of Discretization Error 

Having established convergence of the error indicator, we now examine 
the level of discretization error across the database. Figure 17 shows the 
value of the final error estimate £ in the carpet plot at upper-left. The 

J The initial mesh contains four levels of refinement. 
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error estimates are computed with C = | in Eq. 37 for shock-free cases 
that do not require a limiter (M^ < 0.5), otherwise C — 2 is used. The 
carpet plot shows that £ is small and essentially uniform (£ < 2 counts) 
over most of the subsonic and supersonic cases. Not surprisingly, the 
largest error estimates are obtained in the transonic regime, specifically 
when M 0 c = 0.7, which corroborates with the slow convergence of 77 in 
Fig. 16(b). There is also an increase in £ as M ^ -A 0 that indicates 
some loss of accuracy in the incompressible regime due to the lack of low 
Mach-number preconditioning. This is consistent with results presented 
in [39]. 

Along the perimeter of the carpet plot in Fig. 17, we show convergence 
of Cd for the last three adaptations of representative cases. The error bars 
denote the value of £. In general, the insets show that Cd changes less 
than one count over the final adaptation and £ is decreasing by about a 
factor of two per cycle and brackets the final solution over at least the last 
two cycles. For example, the lower-left inset shows the classic subsonic 
case Moo = 0.3 and a = 0°, where Cd correctly approaches zero as 
H -A 0. The error estimate brackets the expected final solution (Cd = 0 
to plotting accuracy) on the last three meshes. The case at M^ = 0.7 
and a = 8 °, shown at top-right, is the main exception. The value of 
Cd is still changing significantly and the error estimate is just starting 
to tighten. More adaptation cycles (smaller cells) would be required 
to reduce the error estimate further. Taken together, Figs. 16 and 17 
provide the quantitative evidence required for verification of every case 
in the database. 


Meshing Requirements 

Figure 18 shows the number of cells required to achieve the error levels 
presented in Fig. 17, as well as snapshots of the final near-body meshes for 
selected cases. The carpet plot shows that the number of cells varies by 
about a factor of four over the database and peaks when 0.9 < M^ <1.3. 
To give an indication of the wall-clock time needed to construct this 
database, each case requires, on average, about four minutes on one core 
of a laptop computer k . 

The largest mesh is obtained when M ^ = 1.1 and a = 8 ° and contains 
39,500 cells. The top-left inset shows that the high cell count is primarily 
due to a detached bow shock far upstream from the airfoil, which is 
indicated by the cluster of small cells along the inset’s left edge. In 
contrast, Fig. 18 shows that cases with M ^ > 2 require many fewer 
cells. Recall that the refinement pattern is driven by the inner product of 
the adjoint interpolation error (t/^tq — ^tl) with the flow residual error, 
Eq. 36. Upstream of the bow shock, the flow residuals are zero because 
the spatial discretization preserves uniform freestream. Concurrently, 

k MacBook Air (2013) with a 1.7 GHz Intel Core i7 processor and 8GB of memory. 
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Figure 18. Number of cells for all freestream conditions with samples of 
final adapted meshes (near-body views) 


the adjoint solution vanishes downstream of the limiting characteristic 
intersecting the airfoil’s trailing edge. This is similar to Fig. 7(c) of 
the supersonic vortex problem, where the adjoint solution is zero when 
outside the functional’s zone of dependence. Consequently, \tj\h vanishes 
upstream of the bow shock and outside the limiting characteristic, as 
indicated in the insets on the right side of Fig. 18 for cases at M ^ = 1.3 
and Mqo — 3. This confines the refinement pattern to a “diamond” that 
contracts in the cross-flow direction with increasing Mach number. As a 
result, the smallest mesh, containing only 8860 cells, is obtained at the 
highest Mach number, specifically at M^ = 3 and a = 4°. 

As Moo — >* 0, the carpet plot of Fig. 18 shows a moderate increase 
in the number of cells that is consistent with the increase in £ shown 
in Fig. 17. The inset at lower-left shows a typical mesh in the subsonic 
regime. At a given angle of attack, the meshes are self-similar in the 
sense that the refinement regions simply enlarge with decreasing Mach 
number to compensate for the weakening pressure variations. 

Overall, Figs. 17 and 18 offer a convincing demonstration of the 
benefits of adaptive mesh refinement and automatic error control for 
simulation verification. Even in this simple 2D example, constructing 
a single mesh that accommodates all the flow regimes, from smooth 
subsonic flow to the disparate shock systems of transonic and supersonic 
flow, is a daunting task that requires orders of magnitude more cells. 
Moreover, uniform refinement of such a mesh to demonstrate verification 
is not practical. 
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Figure 19. Delta-wing-body model with a partial view of the wind-tunnel 
sting 


7.2 Pressure Signature of a Body in Supersonic Flow 

This example demonstrates verification of a simulation that predicts the 
near-field pressure signature of a lifting configuration in supersonic flow. 
Such simulations are used to determine loudness of sonic booms [69, 70]. 
The key to credible noise analysis is the accurate prediction of the weak 
shocks and expansions, and their interactions, that are generated by low 
sonic-boom vehicles. These pressure disturbances are highly susceptible 
to attenuation by discretization error and consequently, these simulations 
hinge on the mesh [71, 38, 23, 72, 43, 73]. 

The simulation involves a popular wing-body model, originally iden- 
tified as “Model 4” in the experiments of Hunton et al [74] and used 
subsequently in many studies, for example [70, 38, 72, 75, 76]. Figure 19 
shows the details of the geometry, including an axisymmetric fuselage, a 
delta wing with diamond-airfoil cross sections and an approximate wind- 
tunnel sting from [38] that forms a step junction at the base and extends 
downstream. The surface tessellation contains 1.3 million triangles. 

Our goal is to predict the pressure signature on a line located 3.6 
model lengths below the model’s centerline and parallel to freestream. 
The output functional is given by Eq. 21, where we set n = 2 to emphasize 
the importance of accurately capturing the peaks of the signature. The 
freestream Mach number is 1.68 and the angle of attack is set to match 
a Cl of 0.15. This corresponds to one of the experiments performed 
in [74]. Monotonicity of the solution is maintained through use of the 
Barth- Jespersen limiter [63]. 

The left frame of Fig. 20 shows the setup for a half-body simulation. 
The model and the sensor (solid horizontal line) are shown on the symme- 
try plane of the near-field region of the initial mesh. The mesh is rotated 
to approximately align the cells with the freestream Mach-wave angle 
and the cells are stretched, about 2:1:2, in the direction of the wave 
angle and spanwise to improve computational efficiency. The initial mesh 
contains 12 x 12 x 8 cells, and the geometry and sensor are intersected by 
just five cells. In other words, except for the rotation and cell stretching, 
the mesh is not biased to anticipate features of the solution and should 
be far from the asymptotic region of the functional. The refinement 
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Initial Mesh: Mesh after 12 Adaptations: 


Figure 20. Near-field view of the initial mesh (left), mesh after 12 adapta- 
tions (middle) and isobars (right) with inset frame showing isobars close 
to the model. Pressure signature example (symmetry plane, isobar range 
0.65-0.78 in 0.005 increments, = 1.68) 


threshold for each adaptation cycle is set to the mean error value if the 
error indicator is decreasing, otherwise it is set to twice the mean value 
to limit over-refinement. The maximum number of cells, N max , is set to 
20 million cells. 

The middle and right frames of Fig. 20 show the mesh and flow so- 
lution after 12 adaptations. Before addressing verification, we briefly 
discuss the salient features of the refinement pattern. As expected, the 
refinement follows the shocks and expansions between the model and the 
sensor; however, note the refinement above the model and below the sen- 
sor. Similar to the vortex problem of Sec. 6, we visualize the functional’s 
zone of dependence through the adjoint solution. Figure 21 shows the 
absolute value of the first adjoint variable. The grayscale map is tuned to 
indicate where point-source perturbations of the mass-conservation equa- 
tion influence the pressure signature. In addition to the region between 
the body and the sensor, the plot confirms the dependence on regions 
above the model and below the sensor. This is consistent with the ex- 
pected propagation of characteristics in this three-dimensional flowfield. 
Moreover, no refinement occurs upstream of the leading shock because 
the adjoint solution and the flow residual are zero 1 , which makes \t]\h van- 
ish. There is also little refinement in the wake. This happens because the 
adjoint vanishes downstream of the sensor, since residual perturbations 
past the sensor in supersonic flow cannot affect the signature. Lastly, the 
sensitivity of the error indicator to changes in adjoint curvature triggers 
the refinement near the sensor (middle frame of Fig. 20 and Fig. 21) 
to accommodate both the flow and adjoint solutions on the same mesh. 
This is similar to the discussion around Fig. 13 in Sec. 6.2. 

lr The adjoint is zero due to the quadratic form of J. 
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Figure 21. Contours of density- adjoint (l^ol? white 0, black 0 . 012 ) on 
symmetry plane 

To demonstrate verification, Fig. 22 presents convergence of the sig- 
nature, the functional and its error estimates. Figure 22(a) shows the 
classical, qualitative evidence of mesh convergence by examining changes 
in the signature as the mesh is refined. The main features of the signa- 
ture are captured within ten adaptations and the signal is converged to 
plotting accuracy in twelve adaptations. We also validate the signature 
using experimental data (* symbols) from [74]. The agreement between 
the experiment and the simulation is excellent. 

Figures 22(b)-22(d) show quantitative evidence of convergence as the 
mesh is refined. Figure 22(b) shows convergence of the functional Jh 
(solid circles) and the corrected functional J c (open squares) . Changes in 
these quantities begin to taper once the mesh reaches 200,000 cells, where 
the corrected functional begins to predict the value of the functional on 
the next mesh with increasing accuracy. Figure 22(c) shows convergence 
of the various error estimates. In general, the error increases over the first 
five cycles and thereafter decreases. As in the airfoil database example 
(recall Fig. 16), the increasing error indicates new features emerging in 
the flowfield due to the coarse initial mesh. 

The top line (solid circles) in Fig. 22(c) shows convergence of the 
error-indicator bound 77 . The second line (squares) shows convergence 
of the error estimate £ of Eq. 37 with C — 2 . This assumes Jh is 
0(h), which is conservative and expected because both the flow and 
geometry are dominated by non-smooth features. The error-indicator 
bound is larger than the error estimate, but both show similar convergence 
behavior. The third line (solid triangles) monitors the magnitude of the 
functional update (A J = \ J l H — for each cycle of the adaptation. 

Once the mesh reaches 70,000 cells, the update begins to converge and 
becomes smaller than the error estimate. The fourth line (x symbols) 
shows the magnitude of the remaining error term, Eqs. 36 and 38 without 


34 




(a) Pressure signatures, experimental data 
from [74] 


(b) Convergence of functional (Jh) and 
corrected functional ( J c ) 




(c) Convergence of error estimates (d) Functional convergence with error bars 

showing estimate of discretization er- 
ror 

Figure 22. Mesh convergence for the pressure signature example 


localization. Localization increases the magnitude of the error indicator 
77 by about two orders. This is only slightly larger than what is shown in 
Fig. 10(b) for the vortex problem despite the significantly more complex 
flow. Furthermore, this confirms that the remaining error term in Eq. 32 
is small relative to the adjoint correction. 

Figure 22(d) combines information presented in Figs. 22(b) and 22(c) 
to concisely demonstrate simulation verification. The convergence of 
the functional is shown with error bars that represent the estimated 
discretization error £ (taken from the second line of Fig. 22(c)). The 
largest error occurs on the fourth adaptation and the error estimate 
sharpens considerably by the tenth adaptation (600,000 cells). Starting 
from this mesh, the error estimates reliably bracket the solution obtained 
on the finest mesh (29.4 million cells) indicating asymptotic convergence. 

Figure 23 shows error histograms and the mean error value for the 
initial, fifth and final meshes. The left frame of Fig. 23 shows that most 
of the cells of the initial mesh contribute no error, which is due to the 
uniform starting mesh. The middle frame of Fig. 23 shows the histogram 
after five adaptations (sixth mesh), where the error indicator 77 reaches 
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Figure 23. Histograms showing the distribution of error-indicator values 
for the initial mesh (left), mesh after five adaptations (middle) and final 
mesh (right) for the pressure signature example. 


its maximum value (as shown in Fig. 22(c)). This mesh contains about 
11,000 cells and the histogram shows the beginnings of error equidistri- 
bution, with the mode of the distribution already close to the mean. The 
right frame of Fig. 23 presents the histogram of the final mesh. The 
mode moved to the left, the mean stayed close to the mode and the error 
distribution sharpened considerably. 

In terms of performance, the total simulation (wall-clock) time was 
1 hour 42 minutes on 96 Intel Xeon E5-460L cores [77]. The flow and 
adjoint solutions on the final mesh (29.4 million cells) required about 
20 minutes each. In an engineering setting, Fig. 22(d) indicates that it 
would be sensible to terminate the simulation after twelve adaptations 
(omit the last two solutions). Convergence of the functional and error 
becomes predictable by the tenth adaptation and the level of error is quite 
small after twelve adaptations (4.5 million cells). Moreover, Fig. 22(a) 
shows that after twelve adaptations there are essentially no changes in 
the signature. If the simulation is stopped after twelve adaptations, the 
wall-clock time is only 18 minutes, where the flow and adjoint solutions 
require about 4 minutes on the final mesh. Additional speedup is possible 
by omitting the error analysis on the final mesh, which would reduce the 
net wall-clock time to about ten minutes on these CPU’s. 

7.3 Launch- Abort Vehicle 

The final example involves the prediction of aerodynamic forces for a 
launch- abort- vehicle prototype. The example explores the limits of Pro- 
cedure 1 when applied to problems that contain regions of separated, 
recirculating flow. Referring to Figs. 24 and 25, the vehicle consists of 
a crew module with a tower-mounted abort system. During a launch 
emergency, the four Abort Motors (AMs) ignite to pull the crew mod- 
ule safely away from the rocket stack. Throughout the abort, stability 
and control of the vehicle is maintained by differential thrust from eight 
Attitude Control Motors (ACMs) near the nose. For this example, we 
consider a high-altitude abort case studied previously in [78] at = 4 
and a — 20° that involves strong ACM jets. 
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Heatshield 



Figure 24. Launch- abort vehicle surface geometry 



Figure 25. Close-up view of the nose-cone (shown with triangulation) 
and the ACM nozzles. The thrust settings are indicated in the cutaway 
view of the right frame. 


A close-up view of the 
nose with the ACM nozzles 
is shown on the left side of 
Fig. 25. The jet boundary 
conditions are applied at the 
throat faces of the nozzles, 
which are recessed inside the 
tower, as shown in the mid- 
dle frame of Fig. 25. The 
jet boundary conditions match 
the exit momentum, pressure and thrust of the ACM performance 
data [78]. The right side of Fig. 25 shows the thrust setting for each 
nozzle. The largest thrust is generated by the bottom (south facing) 
nozzle and five of the eight nozzles are active. 

The surface discretization of the vehicle contains roughly 380,000 
triangles and is obtained directly from a CAD model. Figure 25 shows 
an example of the triangulation for the nose-cone component. Figure 26 
shows the initial mesh containing 12,880 cells, with only 580 cut cells. 
As in the previous examples, there is no bias to anticipate features of the 
flow solution and multiple ACM nozzles are contained within a single 
cut cell. The refinement threshold for each adaptation cycle is set to the 
mean error value if the error indicator is decreasing, otherwise it is set to 
twice the mean value (same as the pressure signature example 7.2). The 
output of interest is a weighted sum of normal and axial force coefficients 

j = C N + 0 AC a (42) 



Figure 26. Near-body view of initial 
mesh on symmetry plane 
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Figure 27. Near-body views of the mesh after 12 adaptations = 4 
and a — 20°) 


and N max is set to 50 million cells. 

The simulation uses the minmod limiter and runs for 13 adapta- 
tions. For the last two adaptations the mesh growth is fixed at 2.5 and 
3. Figure 27 shows the near-body mesh after 12 adaptations. The mesh 
contains 16.9 million cells with the finest cells located near the nose, in 
the windward region near the centerline and at the heatshield shoulder. 
The refinement of the wake and above the vehicle is moderate, as shown 
in the inset of Fig. 27. Figures 28 and 29 show the solution after 12 
adaptations. The Mach number contours reveal the interaction of several 
shocks and expansions. The left frame of Fig. 29 shows a close-up of the 
nose bow-shock and the south (downward pointing) ACM jet on symme- 
try plane. The bow-shock impinges on the ACM jet-shock and creates 
a channel of supersonic flow that decelerates through a series of shocks 
to subsonic flow just upstream of the ACM nozzle, as shown in the left 
inset of Fig. 29. This is a well-known shock-shock interaction pattern 
that Edney [79] classified as Type IV interference. The right frame of 
Fig. 29 shows the front view of the solution on a plane just behind the 
ACM nozzles. Note the distortion of the bow shock due to the ACM jets. 

The salient feature of the flow is the interaction of the ACM jets with 
the abort motor and crew module surfaces. The jet paths cannot be 
determined in advance of the computation, yet they have a significant 
influence on forces, moments and heat transfer. In particular, the interac- 
tion of the jets with the surrounding flow can adversely affect the ACM’s 
authority over pitching moment [78]. Figures 27 and 28 show that the 
finest cells track the bow and jet shocks away from the body and the 
first transverse cut-plane shows a highly refined mesh for the three main 
jets near the nose. Such fine cells are necessary because discretization 
errors introduced in this upstream region influence the functional over 
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Figure 28. Near-field view of Mach isocontours after 12 adaptations on 
the symmetry plane. Contours are compressed using \[M to improve jet 
visualization. The surface of the vehicle is shaded by pressure coefficient 
(Moo = 4 and a = 20°). 



Figure 29. Close-up views near the nose showing Mach isocontours: side- 
view on symmetry plane (left) and front-view with a cutting plane just 
behind ACM nozzles (right). Inset on left shows details of the ACM 
shock interaction (M^ = 4 and a = 20°). 
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the rest of the vehicle. Further along the vehicle, the second cut-plane 
of Fig. 27 shows the adaptation pattern still tracking the jet paths, but 
with less refinement. The final cut-plane through the heatshield shoulder 
is nearly symmetric, indicating that at this downstream station the jets 
have negligible influence on the functional. 

To further investigate the refinement of the jets, we use iso-surfaces 
of stagnation enthalpy to visualize the jet paths around the vehicle in 
Fig. 30. The jet surfaces are shaded by Mach number. The bifurcation 
of the main jet is striking. This is a consequence of the over-pressure 
at the ACM nozzle exit caused by the Type IV interference, followed 
by downstream interactions of the jet with the supersonic flow. The 
bifurcated jet contacts the sides of the crew module as it is swept upward 
by the high angle-of-attack flow. The “Below” view of Fig. 30 shows that 
the shocks from the abort motor nozzles alter the jet paths away from 
the crew module and thereby alter the heating experienced by the crew 
module. Lastly, as indicated previously in Fig. 27 and shown in Fig. 30, 
the refinement of the jets wanes as the flow approaches the heatshield 
shoulder. While in reality the jets persist well into the wake of the vehicle, 
the error analysis truncates their refinement once they pass beyond the 
heatshield, where they can no longer impact the functional. This is the 
same behavior as we observed in the supersonic regime of the airfoil 
database example, but here in a much more complex flow. 

Figure 31 summarizes the convergence of the functional, the error 
estimates and the aerodynamic forces. Figure 31(a) shows the value of 
the functional on each mesh along with error bars representing the level of 
discretization error £ via Eq. 37 with C — 2. As in the pressure signature 
example, we assume that Jh is 0{h) because the flow is dominated by 
shocks and contact discontinuities. The large error estimates in the 
second and third adaptation are primarily associated with refinement of 
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(a) Functional convergence with error bars showing 
estimate of discretization error 




(b) Convergence of error estimates (c) Convergence of normal and axial 

force coefficients 

Figure 31. Mesh convergence for the launch-abort vehicle example 


the ACM jets. The inset in Fig. 31(a) shows that over the last three 
adaptations the error bars bracket the functional and the changes in the 
functional are less than 1%. Figure 31(b) shows convergence of the error- 
indicator bound 77 , the error estimate £ and magnitude of the functional 
update (A J = | J l H — J^ - 1 1). While convergence is not as convincing as 
in the previous examples, all three error measures are decreasing once 
the mesh reaches about 2 million cells. 

Despite the complicated flowfield, especially the recirculations behind 
the heatshield and the inactive abort motors, Procedure 1 remains effec- 
tive in addressing critical regions of the flow that influence the outputs. 
Figure 31(c) shows that there are virtually no changes in the axial and 
normal force coefficients over the last two adaptations. More broadly, 
this example characterizes the performance of Procedure 1 in verification 
of difficult engineering simulations, and demonstrates the benefits of af- 
fordable and automatic error control in providing insight into complex 
flowfields. 
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8 Conclusions 


A procedure for automatic control of discretization error in steady simu- 
lations of inviscid flow has been presented. The procedure involves error 
estimates based on the method of adjoint- weighted residuals in conjunc- 
tion with incremental mesh enrichment based on a Cartesian cut-cell 
method. We draw the following conclusions from our results: 

• The code verification example demonstrates that both the func- 
tional Jh and corrected functional J c are 0(/i 2 ), which is the ex- 
pected order of accuracy. The error-indicator bound 77 converges 
at about the same rate due to the discontinuous adjoint field and 
localization. This is sufficient for the discretization error estimate 
£ to be 0{h 2 ) with effectivity close to one. 

• The simulation verification examples show that once the mesh is 
sufficiently fine, £ reliably brackets Jh and decreases by about a 
factor of two per adaptation cycle. 

• The error map \t)\h reliably identifies critical regions of the mesh. 
Consequently, the procedure reliably handles extremely coarse ini- 
tial meshes and generates a converging sequence of affordable meshes 
that accurately predict the output. 

• The procedure offers a practical alternative to the manual gener- 
ation of simulation-specific meshes that encompass the knowledge 
and experience of specialists and instead automatically delivers 
meshes and error estimates that provide significant insight into the 
simulation. 

• Since a mesh refinement study is intrinsic to every simulation, the 
procedure automatically provides convergence histories of outputs 
of interest and error estimates. This makes the task of obtaining 
verified simulations straightforward and essentially automatic. 

• Returning to the opening quote of the Introduction; the results 
show that the procedure attains a sufficient level of reliability and 
robustness for routine use in simulation verification. 

There are many areas for future work. For example, improving the ac- 
curacy of the embedded adjoint representation through use of affordable 
approximate solutions could improve the sharpness of the error estimate. 
Incorporating the procedure within simulation-based design could sig- 
nificantly reduce cost of optimization problems. More broadly, since 
simulations in practice frequently involve some degree of unsteady flow 
(as in the launch abort example of Sec. 7.3), the procedure should be 
extended to transition smoothly from steady to unsteady flow, especially 
for applications where the time variation in the outputs is relatively small. 
Furthermore, the extension of the procedure to handle flows involving 
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multiple components, or species, and turbulence is an important area of 
future work. 
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